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Abstract 

We derive analytical shape derivative formulas of the system mat- 
rix representing electric field integral equation discretized with Raviart- 
Thomas basis functions. The arising integrals are easy to compute with 
similar methods as the entries of the original system matrix. The results 
are compared to derivatives computed with automatic differentiation tech- 
nique and finite differences, and are found to be in excellent agreement. 

1 Introduction 

The adjoint variable methods based on shape gradients have recently gained 
significant attention in the shape optimization of microwave devices [fJl!15| . The 
main reason for such an interest is that the adjoint variable approach eases sig- 
nificantly the computational burden of the optimization process: computing the 
gradient of the objective requires only one solve of the governing state equations, 
and in addition, at most one solve of the adjoint problem. 

The adjoint variable methods rely on the availability of the derivatives of 
the system matrix with respect to the control parameters. Traditionally these 
derivatives have been computed with finite difference (FD) formulas. The main 
drawback of this approach is that the accuracy greatly depends on the step 
length parameter, which has to be carefully chosen in order to find a suitable 
balance between round off and truncation errors. 

Use of the automatic differentiation (AD) was proposed in [15] to compute 
the derivatives of the system matrix arising from a method of moments dis- 
cretization of the electric field integral equation (EFIE). The method was ap- 
plied to antenna shape optimization problems in [151116] . AD computes the exact 
derivatives of the computer code, but there are some complications involved. 
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AD can be implemented either as a source transformation tool, or using op- 
erator overloading. Source transformation tools are complicated to implement, 
and existing tools are available only for some programming languages. While 
using these tools, one may be forced to restructure the code, and avoid such fea- 
tures of the programming language that are not supported. Tools based on the 
operator overloading are more simple to implement, but not all programming 
languages support operator overloading. Moreover, there is some execution 
overhead related to this approach. In both cases, AD has to be applied also to 
all subroutines that are called from the code, which may present difficulties if 
one wishes to use any external libraries. 

The electric field integral equation is a widely used method to compute scat- 
tering of time harmonic electromagnetic field from perfectly electrically conduct- 
ing (PEC) bodies and surfaces [3,8. 10J. In this paper we derive simple analytical 
formulas for calculating the system matrix derivatives of the discretized EFIE 
system. The arising integrals are such that they are easy to compute using 
methods similar to those used to calculate the elements of the original system 
matrix. The results are compared against derivatives calculated using automatic 
differentiation and difference formulas. 

We employ the lowest order Raviart-Thomas (RT) basis functions [13] in 
the discretization. They are almost identical to the Rao-Wilton-Glisson (RWG) 
basis functions [T2], which are the usual choice to discretize the EFIE, except 
that usually the RWG functions are scaled with the length of the edge which 
they are related to. 

The shape derivatives of the electromagnetic field solution together with the 
far field pattern has been characterized and analyzed in [4][5l[TT]. Contrary to 
these works, we consider differentiation of the discretized system, and utilize 
the adjoint variable method. 

2 Preliminaries 

Let S be a surface with or without boundary in R 3 and T its triangulation. It 
can be a boundary of a some sufficiently regular open domain in R 3 , in which 
case it models the scattering object. If it is not a boundary of any set, but it 
otherwise exhibits sufficient regularity, it acts as a model for a thin perfectly 
conducting screen. 

In what follows, RT(T) denotes the space spanned by the lowest order 
Raviart-Thomas (RT) basis functions on T, V s • denotes the surface divergence, 
and V s denotes the surface gradient. 

In the EFIE the unknown function to be solved is, roughly speaking, the 
equivalent surface current J = nxH, where n is the surface unit normal and 
H is the total magnetic field. The source term in the equation is the incident 
electric field denoted by E p . 

The current J satisfies the Rumsey reaction principle [Tl 121181m] : 
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Find J e H x 1/2 (y s -;S) s.t. 

— / V s • u(x) / g k (x - y)V s ■ J(y)dydx - iu/i, / u(x) • / g k (x - y)J(y)dydx 
we Js Js Js Js 

= I u(x) ■ B p (x)dx Vu 6 H~ 1/2 (V a -;S). (1) 
Js 

Here is the fundamental solution of the Hehnholtz equation: 

e ik\x\ 

9k(x) = -— w (2) 

and V s - is the surface divergence. We denote the Euclidean norm of x G M 3 by 
\x\. 

For detailed treatment of the space H x ' (V S -;S) and the above equation 
we refer to [2]. We note that the integrals in ((!]) should be interpreted as a 
certain duality pairing, but at the discrete level they are proper integrals. 

After discretizing the equation ([T]) we arrive to the following finite dimen- 
sional variational problem: 

Find Jh 6 RT such that 

a(v, J h ) = / v(x) ■ E p (x)dx, Vv e RT(T), (3) 
Js 

where the bilinear form a is given by 

a(v,u)=— / V s -v(x) / g k (x - y)V s • u(y)dydx- 
we Js Js 



iuj/j, / v(x) • / gk(x - y)u(y)dydx. (4) 

We denote the corresponding linear system of equations of ([3]) by 

Ax = b. (5) 

Let T p ,Tg £ 7", u, v £ RT(7~). and suppu and suppv intersect T p and T q , 
resp. In order to compute the entries of A one needs to compute integrals of 
type 

7i(u,v;T p ,T g ) = / v(x) ■ / g k {x - y)\i{y)dydx, (6) 

I 2 (u,v;T p ,T q )= [ V,-v(a)/ g fc (x - y)V s • n{y)dydx. (7) 

We denote the usual interpolating nodal piecewise first order polynomials 
on 5 associated with T by (A m )^ =1 , where N is the number of vertices in 1~. 
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2.1 Adjoint variable methods for the EFIE system 



Let a be the vector of design variables, and J7"(a,x(a)) be a real valued objec- 
tive function. In general, real valued functions are not differentiable with respect 
to complex arguments in the conventional complex analytic sense. Therefore we 
differentiate J separately with respect to the real and imaginary parts of the 
variables x: 

dj_ _ dj_ dj_mx + dj_d^ 

dak dat (?3?x dak <9Sx dak 

Here dj /dak reflects the explicit dependence of J on the design. 

The convention V x J := + iVQ xt 7 will be used to simplify notation. 

It can be easily checked that equation (|5]) can now be written as 



dj 

dak 



dj 

da k 



(V X J) H 



<9x 

da k 



By differentiating the state relation Ax = b, we obtain 



<9x 

dak 



dA 

da k 



db 

da k 



(9) 



(10) 



In the so called direct differentiation approach, the derivatives dx/dak are solved 
from (|10p . and is used to compute the gradient of the objective function. 
However, the right hand side of (1101) is different for each design variable, which 
makes this approach rather inefficient. 
By introducing the adjoint problem 



A H -/ = V x J, 

and using (fT0|) . equation (j9]) can be written as 



dj 
da k 



dj 
da k 
dj 
da k 
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\dak dak 

{i h _ dA^ 

\da k da k 



(11) 

(12) 
(13) 



In the adjoint variable method, equation (fTTj) is solved for 7, and equation (jT5|) is 
used to compute the gradient of the objective. Notice that the adjoint problem 
(TTTj) does not depend on the design, and therefore 7 is the same for all design 
variables. In case of some particular objective functions (see e.g. [9]), the adjoint 
vector 7 can be immediately obtained from the solution vector x, and one does 
not have to solve the adjoint problem at all. 

In any case, the derivatives of the right hand side and the system matrix 
with respect to the design variables are required. This paper presents analytical 
formulas, which can be used to compute them efficiently in the case where A is 
a system matrix arising from the EFIE formulation ©. 
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3 The derivative formulas 



We define the shape deforming mapping with the aid of the global nodal shape 
functions A m as follows. 

Let r £ I 3 be a vector and < s < 1. We define the deformation mapping 
FJ™ associated with node m by 

F™(x) =x + S T\ m (x). (14) 

Note that the domain of F™ is the whole of S. We shall drop the superscript 
m from F™ whenever it does not cause confusion. 

Let K be a triangle in M 2 and F a diffeomorphism from K to its image. We 
define the Piola transform Vf associated with F by 

P F u=^FW- (15) 

which maps smooth sections of the tangent bundle of K to those of F{K). Here 
the Jacobian of F by is denoted by F' . We note that the divergence of the Piola 
transformed function u is given by 

v *'^ u= dcT^ (Vs ' u)oF_1 - (16) 

This simplifies the calculations greatly, as we shall see. 

The Raviart-Thomas basis functions can be defined using the Piola transform 
and reference element as follows. Let K be the standard 2-simplex 

{(x,y) EM 2 : < x, < y, andO < x + y < 1} 

and define functions on K by 

ui(z,y) = [x, 1 - y] T , 

U2{x,y) = [x, y] T and (17) 

u 3 (x,y) = [1 -x, y] T . 

Now the restriction of an arbitrary u G RT(T) to the triangle if is a Piola 
transform of some ie. 

u = V Fk Si 

for some i € {1,2,3}, where Fk is the affine mapping which maps K to K . We 
note that the the determinant of F' K is given by detF^- = ±2|if|, where \K\ is 
the area of K and the sign is determined by its orientation. 
The main objects of interest in this work are the derivatives 

^-I m (V Fs u,V Fs v;F s (T p ),F s (T q ))\ s=0 , m = l,2. (18) 



•5 



To start with, we note that 
h {T Fs u, V Fs v; F s (T p ), F s (T,)) = 



■ , r', > l i t 1 // \ 9k{x - y)dydx = 

F s (T q ) J f 3 (t p ) detF^x) deti^(y) 

F»v(z) • F^y)u(y)g k (F s (x) - F s (y))dydx. (19) 



and 



/ 2 (Pf s u, 7> Fs v; F s (T p ),F s (T q )) = 

(V s -v)oF- 1 (x) (V.-u)oi?-i( v ) 



Qi-fx — y)dydx = 
Fs (T q )JF s (T p) detFKx) detFj(y) yu yj 

/ / V s -v(x)-V s -u(y) 5fc (F s (x)-F s (y))rf^a;. (20) 

For conciseness, we denote = / m ("Pf s u, 'Pf, v ; Fs(T p ), F s (T q )). 

Looking at (fl9|) and (|20| we find that we only need to compute derivatives 
of gk and F' a wrt. s. 

Lemma 1. Let F™ oe given by Eq. (fT4|) . If /io/ds i/iai 

£1*™ -™u = ^yti^Mz^m (21) 



and 



^F s '(x) T F s '(y)| s=0 = (T\7 s X m (x)) T +rV s X m (y)- (22) 



Proof. Just by calculating 

- Fr(y)| 2 = \x - y + sr (X m (x) X m (y)) | 2 

= A + 2s{x - y) ■ r(X m (x) - X m (y)) + Bs 
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where A and B are constants wrt. s. By using chain rule, the first assertion 
follows. 

For the second formula, we note that 

F s '(z) =I+STV s X m (x) (23) 
and the assertion is obvious. □ 
Using the above results we have the shape derivative of the kernel gk as 

(24) 
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Furthermore, this can be simplified to 
d 

—g k (F s (x) - F a (y))\ a=0 = r • (Vg k ){x - y)(X m (x) - X m (yj). 
Thus, we immediately get the formula 



= ( ( ■ v(x)V s • u(j/)r • {Vgk){x -y) (X m (x) - X m (y)) dydx. 



(25) 

The derivative Jj/f is only slightly more complicated; it is obtained by an 
application of the Leibniz's rule: 

Ts* 1 =11 (§~s Ax) ' F s( x f F s(vMy)) 9k{x - y)+ 

v(x) • n(y)^-g k {F s {x) - F s (y))dydx. (26) 

Thus we obtain 



=11 v ^ ' {( T ^sX m (x)) T + TV s X m (y)) ■ u(y)g k (x - y)+ 
v(x) • u(y)r • {Vg k ){x - y) (X m (x) - X m (y)) dydx. 



(27) 



4 Numerical verifications 



We compare the results given by the formulas (1271) and (|25| to the values of 
derivative computed with difference formulas and automatic differentiation ap- 
plied to numerical code which computes at s = 0. 

To start with, we note that the integrals (|2"T)) and ([23)1 are singular when 
the closures of T p and T q intersect. Thus we manipulate the expressions in such 
a way that we can directly apply singularity subtraction methods to compute 
them. 

The singularity subtraction method is designed to calculate singular integrals 
of the form 

' v(x) J K(x - y)u(y)dydx, (28) 

where K is singular at 0. The idea is to express K as a sum of singular K s and a 
function Kb bounded at 0, in such a way that J K s (x— y)v(x)dx can be evaluated 
analytically, and J Kb(x—y)v(x)dx integrates easily with numerical quadratures. 
In this paper, we use the same quadrature to integrate K = K s + Kb wrt. y 
and the inner integral with kernel Kb wrt. x. Note that J Kb(x — y)v(x)dx is 
integrated analytically. 
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Looking at equations (|25p and (|2"Tj) we observe that we only need to compute 
integrals of type 



V s • v(x) / t ■ V x g k (x - y) (X m (x) - \ m (y)) V s • u(y)dydx 
v(x) ■ J t ■ V x g k (x - y) (A m (x) - A m (t/)) u(y)dydx (29) 
v(x) • T(x,y)u(y)g k (x - y)dydx 



Here T(x,y) = (rV s A m (a;)) T + rV s A m (y). Looking at these, we find that 
they can be calculated using singularity subtraction techniques discussed e.g. 

In the following we shall compare numerical values of 

§- s I s m \ s=0 , m=l,2 (30) 

calculated with automatic differentiation applied to the singularity subtraction 
technique, analytical formulas introduced here, and forward difference formula. 
We inspect these values with four different triangle configurations that occur in 
computations: the triangles T p and T q do not touch, they share a vertex, they 
share an edge, and the triangles are equal. These four cases are presented in 
Figure [TJ 

The automatic differentiation algorithm generates code which computes the 
derivatives of (|30|) . The non-differentiated code calculates the system matrix 
contributions of two triangles with singularity subtraction technique by remov- 
ing the two most singular terms from the kernel of the EFIE integral. To get 
directly comparable results, we also subtract two terms when computing the 
derivatives using analytical formulas. For details and an example where the AD 
technique is employed in EFIE calculations we refer to [15J. 

We compute the derivatives of the local system matrix contributions 

o^(u,v) =^- £ h{V Fs u,V Fa v;F s (T p ),F s {T q ))- 

iuiih{V Fa u, V Fs v; F s (T p ), F s (T q )), (31) 

where u and v vary over three local RT basis functions on T p and T q respectively, 
and investigate the effect of changing the number of integration points n. The 
quadrature rules were adapted from [B]. 
In Tables [TaHldl we inspect the sum 

3 E^<( U ' V )' (32) 

U. V 

computed with the AD method, analytical formulas and the FD method given 

by 

f'(s)*(f(s + h)-f(s))/h. (33) 
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The step length parameter was chosen to be h = 10~ 8 . The reason for presenting 
only the imaginary part is that it is the most involved integral to compute. 

In Tables l2aH2dl we have the difference in the Frobenius norm of the 3x3 
matrices -§^o- s pq relative to the values computed using the present methods. 

We see that the proposed analytical formulas and the automatic differenti- 
ation approach produce same results up to numerical precision. Relative differ- 
ence of the derivatives obtained with the FD method is of the order 1CP 7 . This 
is expected, since the forward finite difference formula has a truncation error of 
the order 0(h). 

5 Concluding remarks 

We have derived compact analytical formulas to calculate the shape derivatives 
of the EFIE system matrix arising from RT discretization. The derivatives are 
easy to implement, since the emerging integrals can be calculated using similar 
algorithms as the original EFIE system matrix. In this case we used singularity 
subtraction technique, but the use of the singularity cancellation methods |17| 
could be possible as well. 

We have shown that the shape derivatives are the same, up to numerical 
precision, as the ones computed using automatic differentiation (AD), if one 
uses the same techniques to evaluate the arising integrals. On the other hand, 
AD is known to produce the exact derivatives of the given computer realiza- 
tion. Such derivatives are perfectly consistent with the objective function val- 
ues, which is very convenient from the optimization perspective. The accuracy 
of the derivatives is naturally dictated by the utilized numerical methods, and 
depends heavily for example on the number of points used for the numerical 
integration. 

The analytical formulas offer several benefits over the automatic differen- 
tiation approach. Existing AD tools are available only for some programming 
languages, and while using these tools, one may be forced to restructure the code 
and avoid some features of the programming language. Moreover, AD has to be 
applied also to all subroutines that are called from the code, which may present 
difficulties if one wishes to use any external libraries. Finally, implementing the 
analytical formula is likely to lead to better computational performance, as one 
avoids execution overhead related to some AD techniques, and the programmer 
is free to optimize the code. 

The FD derivatives are approximate by nature, and their accuracy greatly 
depends on the step length parameter. From the numerical results we can con- 
clude that the derivatives computed with the proposed formulas are in agreement 
with the finite difference derivatives. Use of the analytical formulas is therefore 
preferred, because the need to choose a step length parameter is avoided. 

By having an analytical formula at hand, one is free to use any suitable 
method for its numerical evaluation. Thus it is easy to make a trade-off be- 
tween accuracy and computational efficiency. Moreover, the analytical formulas 
provide a basis for understanding the behaviour of the derivatives. 
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(a) near (b) point (c) edge (d) same 



Figure 1: Three different configurations of triangles T p and T q triangles. The 
triangle T q remains still in each of the cases. 



Table 1: The imaginary part of I(u, v) defined by the equation (|32[) 

(a) near 



n 


Analytical 


AD 


FD 


6 


2 


98169473e+01 


2 


98169473e+01 


2 


98169583e+01 


16 


2 


98097099e+01 


2 


98097099e+01 


2 


98097110e+01 


61 


2 


98096873e+01 


2 


98096873e+01 


2 


98096885e+01 


85 


2 


98096873e+01 


2 


98096873e+01 


2 


98096876e+01 



(b) point 



n 


Analytical 


AD 


FD 


6 


-6 


13463371e+00 


-6 


. 13463371e+00 


-6 


13463194e+00 


16 


-6 


11420678e+00 


-6 


. 11420678e+00 


-6 


11420794e+00 


61 


-6 


11421020e+00 


-6 


. 11421020e+00 


-6 


11421161e+00 


85 


-6 


11416794e+00 


-6 


. 11416794e+00 


-6 


11417224e+00 










(c) edge 






n 




Analytical 




AD 




FD 


6 


-7 


84966165e+01 


-7 


.84966165e+01 


-7 


84966142e+01 


16 


-7 


67762767e+01 


-7 


.67762767e+01 


-7 


67762887e+01 


61 


-7 


68970360e+01 


-7 


.68970360e+01 


-7 


68970451e+01 


85 


-7 


68875833e+01 


-7 


.68875833e+01 


-7 


68875768e+01 










'd) same 






ii 




Analytical 




AD 




FD 


6 


-5 


98000685e+01 


-5 


.98000685e+01 


-5 


98000312e+01 


16 


-5 


86918000e+01 


-5 


.86918000e+01 


-5 


86917977e+01 


61 


-5 


87792043e+01 


-5 


.87792043e+01 


-5 


87792522e+01 


85 


-5 


87738744e+01 


-5 


.87738744e+01 


-5 


87738958e+01 
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Table 2: Relative difference to analytical derivative 

(a) near (b) point 



n 


AD 




FD 




n 


AD 


FD 


6 


1 


10e-15 


2 


14e- 


07 


6 


1 


07e- 


15 


1 


30e- 


07 


16 


1 


lle-15 


5 


27e- 


08 


16 


7 


25e- 


16 


9 


42e- 


08 


61 


1 


85e-15 


3 


02e- 


08 


61 


2 


05e- 


15 


4 


73e- 


08 


85 


2 


52e-15 


2 


74e- 


08 


85 


1 


24e- 


15 


1 


57e- 


07 






(c) edge 










(d) 


same 






n 




AD 




FD 




n 




AD 






FD 




6 


5 


53e-16 


3 


51e- 


08 


6 


6 


13e- 


15 


4 


17e- 


07 


16 


3 


21e-15 


1 


16e- 


07 


16 


1 


84e- 


15 


9 


87e- 


08 


61 


7 


59e-15 


9 


78e- 


08 


61 


3 


06e- 


15 


5 


39e- 


07 


85 


7 


82e-15 


6 


34e- 


08 


85 


1 


08e- 


15 


2 


44e- 


07 



References 

[1] A. Buffa, M. Costabel, and C. Schwab. Boundary element methods for 
Maxwell's equations on non-smooth domains. Numerische Mathematik, 
92:679-710, 2002. 

[2] Annalisa Buffa and Raff Hiptmair. Galerkin boundary element methods for 
electromagnetic scattering. In Topics in computational wave propagation, 
volume 31 of Led. Notes Comput. Sci. Eng., pages 83-124. Springer, Berlin, 
2003. 

[3] Snorre H. Christiansen. Discrete Fredholm properties and convergence es- 
timates for the electric field integral equation. Math. Comput., 73:143-167, 
January 2004. 

[4] Martin Costabel and Fred'erique Le Louer. Shape derivatives of boundary 
integral operators in electromagnetic scattering. arXiv, 2010. 

[5] Martin Costabel and Fred'erique Le Louer. Shape derivatives of boundary 
integral operators in electromagnetic scattering. Part II: Application to 
scattering by a homogeneous dielectric obstacle. arXiv, 2011. 

[6] D. A. Dunavant. High degree efficient symmetrical Gaussian quadrature 
rules for the triangle. International Journal for Numerical Methods in En- 
gineering, 21(6):1129-1148, 1985. 

[7] Seppo Jarvenpaa, Matti Taskinen, and Pasi Yla-Oijala. Singularity subtrac- 
tion technique for high-order polynomial vector basis functions on planar 
triangles. IEEE Transactions on Antennas and Propagation, 54(l):42-49, 
January 2006. 



11 



[8] B.M. Kolundzija and A.R. Djordjevic. Electromagnetic modeling of com- 
posite metallic and dielectric structures. Artech House Electromagnetic 
Analysis Series. Artech House, 2002. 

[9] N.K. Nikolova, Jiang Zhu, Dongying Li, M.H. Bakr, and J.W. Bandler. 
Sensitivity analysis of network parameters with electromagnetic frequency- 
domain simulators. IEEE Transactions on Microwave Theory and Tech- 
nique, 54(2):670 - 681, February 2006. 

[10] A.J. Poggio and E.K. Miller. Computer techniques for electromagnetics, 
chapter Integral equation solutions of three-dimensional scattering prob- 
lems. Oxford, UK: Pergamon Press, 1973. 

[11] Roland Potthast. Domain derivatives in electromagnetic scattering. Math. 
Methods Appl. Sci., 19(15):1157-1175, 1996. 

[12] Sadasiva M. Rao, Donald R. Wilton, and Allen W. Glisson. Electromag- 
netic scattering by surfaces of arbitrary shape. IEEE Transactions on An- 
tennas and Propagation, 30(3):409-418, May 1982. 

[13] P.-A. Raviart and J. M. Thomas. A mixed finite element method for 2nd 
order elliptic problems. In Mathematical aspects of finite element methods 
(Proc. Conf., Consiglio Naz. delle Ricerche (C.N.R.), Rome, 1975), pages 
292-315. Lecture Notes in Math., Vol. 606. Springer, Berlin, 1977. 

[14] V. H. Rumsey. Reaction Concept in Electromagnetic Theory. Phys. Rev., 
94:1483-1491, June 1954. 

[15] J. I. Toivanen, R. A. E. Makinen, S. Jarvenpaa, P. Yla-Oijala, and J. Ra- 
hola. Electromagnetic sensitivity analysis and shape optimization using 
method of moments and automatic differentiation. IEEE Transactions on 
Antennas Propagation, 57(1):168-175, 2009. 

[16] Jukka I. Toivanen, Raino A. E. Makinen, Jussi Rahola, Seppo Jarvenpaa, 
and Pasi Yla-Oijala. Gradient-based shape optimisation of ultra-wideband 
antennas parameterised using splines. IET Microwaves, Antennas and 
Propagation, 4:1406-1414, 2010. 

[17] D. Wilton, S. Rao, A. Glisson, D. Schaubert, O. Al-Bundak, and C. Butler. 
Potential integrals for uniform and linear source distributions on polygonal 
and polyhedral domains. IEEE Transactions on Antennas and Propagation, 
32(3):276 - 281, March 1984. 



12 



